#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "THnSparse.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"

#include "AliMCEventHandler.h"
#include "AliMCEvent.h"

#include <AliCentrality.h>

#include "AliESDtrackCuts.h"
#include "AliStack.h"

#include "AliAnalysisTaskPhi.h"

//#define PT_MIN 0.3  //set range for Pt
//#define PT_MAX 1.5  //(if this selection is performed in ESD track loop)

#define ETA_MIN -0.8  //set range for pseudorapidity
#define ETA_MAX 0.8   //(if this selection is performed in ESD track loop)


ClassImp(AliAnalysisTaskPhi)

//________________________________________________________________________
AliAnalysisTaskPhi::AliAnalysisTaskPhi(const char *name)
    : AliAnalysisTaskSE(name)
    , fVxMax(0)
    , fVyMax(0)
    , fVzMax(0)
    , fOutput(0)
    , fHistPtPlus(0)
    , fHistPtMinus(0)
    , fHistPt(0)
    , fHistPtPlusAfterPtCut(0)
    , fHistPtMinusAfterPtCut(0)
    , fHistPtAfterPtCut(0)
    , fHistEtaPlus(0)
    , fHistEtaMinus(0)
    , fHistEta(0)
    , fHistEtaPlusAfterPtCut(0)
    , fHistEtaMinusAfterPtCut(0)
    , fHistEtaAfterPtCut(0)
    , fHistEtaPlusAfterCuts(0)
    , fHistEtaMinusAfterCuts(0)
    , fHistEtaAfterCuts(0)
    , fHistPtPlusAfterCuts(0)
    , fHistPtMinusAfterCuts(0)
    , fHistPtAfterCuts(0)
    , ffHistEtaPtPlusAfterCuts(0)
    , ffHistEtaPtMinusAfterCuts(0)
    , ffHistEtaPtAfterCuts(0)
    , fTrackCuts(0)
    , fHistCharge(0)
    //, fHistY(0)
    //, fHistCent(0)
    , fHistVtxX(0)
    , fHistVtxY(0)
    , fHistVtxZ(0)
    , fHistTPCNcls(0)
    , fHistITSNcls(0)
    , ffHistDCAxyz(0)
    , fHistChi2PerClusterTPC(0)
    , fHistNch(0)
    , fHistNchAccepted(0)
    , ffHistEtaTPCNcls(0)
    , ffHistEtaPhiPlus(0)
    , ffHistEtaPhiMinus(0)
    , ffHistEtaPhi(0)
    , ffHistEtaPhiPlus64(0)
    , ffHistEtaPhiMinus64(0)
    , ffHistEtaPhi64(0)
    , fHistsSparse(0)
    , ffHistNpEvent(0)
    , ffHistNmEvent(0)
    //, fHistPtEtaEvent(0)
    , fHistNchEvent(0)
    , fHistNchAcceptedEvent(0)
    , fHistEventCutStats(0)
    , fHistTrackCutStats(0)
{
    // Constructor

    // Define input and output slots here
    // Input slot #0 works with a TChain
    DefineInput(0, TChain::Class());
    // Output slot #0 id reserved by the base class for AOD
    // Output slot #1 writes into a TH1 container
    DefineOutput(1, TList::Class());
}

AliAnalysisTaskPhi::AliAnalysisTaskPhi()
    : AliAnalysisTaskSE("user_task")
    , fVxMax(0)
    , fVyMax(0)
    , fVzMax(0)
    , fOutput(0)
    , fHistPtPlus(0)
    , fHistPtMinus(0)
    , fHistPt(0)
    , fHistPtPlusAfterPtCut(0)
    , fHistPtMinusAfterPtCut(0)
    , fHistPtAfterPtCut(0)
    , fHistEtaPlus(0)
    , fHistEtaMinus(0)
    , fHistEta(0)
    , fHistEtaPlusAfterPtCut(0)
    , fHistEtaMinusAfterPtCut(0)
    , fHistEtaAfterPtCut(0)
    , fHistEtaPlusAfterCuts(0)
    , fHistEtaMinusAfterCuts(0)
    , fHistEtaAfterCuts(0)
    , fHistPtPlusAfterCuts(0)
    , fHistPtMinusAfterCuts(0)
    , fHistPtAfterCuts(0)
    , ffHistEtaPtPlusAfterCuts(0)
    , ffHistEtaPtMinusAfterCuts(0)
    , ffHistEtaPtAfterCuts(0)
    , fTrackCuts(0)
    , fHistCharge(0)
    //, fHistY(0)
    //, fHistCent(0)
    , fHistVtxX(0)
    , fHistVtxY(0)
    , fHistVtxZ(0)
    , fHistTPCNcls(0)
    , fHistITSNcls(0)
    , ffHistDCAxyz(0)
    , fHistChi2PerClusterTPC(0)
    , fHistNch(0)
    , fHistNchAccepted(0)
    , ffHistEtaTPCNcls(0)
    , ffHistEtaPhiPlus(0)
    , ffHistEtaPhiMinus(0)
    , ffHistEtaPhi(0)
    , ffHistEtaPhiPlus64(0)
    , ffHistEtaPhiMinus64(0)
    , ffHistEtaPhi64(0)
    , fHistsSparse(0)
    , ffHistNpEvent(0)
    , ffHistNmEvent(0)
    //, fHistPtEtaEvent(0)
    , fHistNchEvent(0)
    , fHistNchAcceptedEvent(0)
    , fHistEventCutStats(0)
    , fHistTrackCutStats(0)
{
    // Constructor

    // Define input and output slots here
    // Input slot #0 works with a TChain
    DefineInput(0, TChain::Class());
    // Output slot #0 id reserved by the base class for AOD
    // Output slot #1 writes into a TH1 container
    DefineOutput(1, TList::Class());
}

AliAnalysisTaskPhi::~AliAnalysisTaskPhi()
{
    // Destructor
}

//________________________________________________________________________
void AliAnalysisTaskPhi::UserCreateOutputObjects()
{
    // Create histograms
    // Called once
    //
    // For histograms description see .h file

    fOutput = new TList();

    fHistPtPlus = new TH1F("fPtPlus", "+ Particles Vs Pt", 1000, 0, 10 );
    fHistPtPlus->GetXaxis()->SetTitle("GeV");
    fHistPtPlus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtPlus);

    fHistPtMinus = new TH1F("fPtMinus", "- Particles Vs Pt", 1000, 0, 10 );
    fHistPtMinus->GetXaxis()->SetTitle("GeV");
    fHistPtMinus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtMinus);

    fHistPt = new TH1F("fPt", "ch Particles Vs Pt", 1000, 0, 10 );
    fHistPt->GetXaxis()->SetTitle("GeV");
    fHistPt->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPt);

    fHistPtPlusAfterPtCut = new TH1F("fPtPlusAfterPtCut", "+ Particles Vs Pt after Pt cut", 600, 0, 6 );
    fHistPtPlusAfterPtCut->GetXaxis()->SetTitle("GeV");
    fHistPtPlusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtPlusAfterPtCut);

    fHistPtMinusAfterPtCut = new TH1F("fPtMinusAfterPtCut", "- Particles Vs Pt after Pt cut", 600, 0, 6 );
    fHistPtMinusAfterPtCut->GetXaxis()->SetTitle("GeV");
    fHistPtMinusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtMinusAfterPtCut);

    fHistPtAfterPtCut = new TH1F("fPtAfterPtCut", "ch Particles Vs Pt after Pt cut", 600, 0, 6 );
    fHistPtAfterPtCut->GetXaxis()->SetTitle("GeV");
    fHistPtAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtAfterPtCut);

    fHistEtaPlus = new TH1F("fEtaPlus", "+ particles vs Eta", 400, -2, 2 );
    fHistEtaPlus->GetXaxis()->SetTitle("#eta");
    fHistEtaPlus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaPlus);

    fHistEtaMinus = new TH1F("fEtaMinus", "- particles vs Eta", 400, -2, 2 );
    fHistEtaMinus->GetXaxis()->SetTitle("#eta");
    fHistEtaMinus->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaMinus);

    fHistEta = new TH1F("fEta", "ch particles vs Eta", 400, -2, 2 );
    fHistEta->GetXaxis()->SetTitle("#eta");
    fHistEta->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEta);

    fHistEtaPlusAfterPtCut = new TH1F("fEtaPlusAfterPtCut", "+ particles vs Eta after Pt Cut", 400, -2, 2 );
    fHistEtaPlusAfterPtCut->GetXaxis()->SetTitle("#eta");
    fHistEtaPlusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaPlusAfterPtCut);

    fHistEtaMinusAfterPtCut = new TH1F("fEtaMinusAfterPtCut", "- particles vs Eta after Pt Cut", 400, -2, 2 );
    fHistEtaMinusAfterPtCut->GetXaxis()->SetTitle("#eta");
    fHistEtaMinusAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaMinusAfterPtCut);

    fHistEtaAfterPtCut = new TH1F("fEtaAfterPtCut", "ch particles vs Eta after Pt Cut", 400, -2, 2 );
    fHistEtaAfterPtCut->GetXaxis()->SetTitle("#eta");
    fHistEtaAfterPtCut->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaAfterPtCut);

    fHistEtaPlusAfterCuts = new TH1F("fEtaPlusAfterCuts", "+ particles vs Eta after all cuts", 400, -2, 2 );
    fHistEtaPlusAfterCuts->GetXaxis()->SetTitle("#eta");
    fHistEtaPlusAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaPlusAfterCuts);

    fHistEtaMinusAfterCuts = new TH1F("fEtaMinusAfterCuts", "- particles vs Eta after all cuts", 400, -2, 2 );
    fHistEtaMinusAfterCuts->GetXaxis()->SetTitle("#eta");
    fHistEtaMinusAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaMinusAfterCuts);

    fHistEtaAfterCuts = new TH1F("fEtaAfterCuts", "ch particles vs Eta after all cuts", 400, -2, 2 );
    fHistEtaAfterCuts->GetXaxis()->SetTitle("#eta");
    fHistEtaAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistEtaAfterCuts);

    fHistPtPlusAfterCuts = new TH1F("fPtPlusAfterCuts", "+ Particles Vs Pt after all cuts", 600, 0, 6 );
    fHistPtPlusAfterCuts->GetXaxis()->SetTitle("GeV");
    fHistPtPlusAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtPlusAfterCuts);

    fHistPtMinusAfterCuts = new TH1F("fPtMinusAfterCuts", "- Particles Vs Pt after all cuts", 600, 0, 6 );
    fHistPtMinusAfterCuts->GetXaxis()->SetTitle("GeV");
    fHistPtMinusAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtMinusAfterCuts);

    fHistPtAfterCuts = new TH1F("fPtAfterCuts", "ch Particles Vs Pt after all cuts", 600, 0, 6 );
    fHistPtAfterCuts->GetXaxis()->SetTitle("GeV");
    fHistPtAfterCuts->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistPtAfterCuts);

    ffHistEtaPtPlusAfterCuts = new TH2F("ffEtaPtPlusAfterCuts", "+ Particles Vs Eta and Pt after all cuts", 200, -1, 1, 600, 0, 6);
    ffHistEtaPtPlusAfterCuts->GetXaxis()->SetTitle("#eta");
    ffHistEtaPtPlusAfterCuts->GetYaxis()->SetTitle("GeV");
    fOutput->Add(ffHistEtaPtPlusAfterCuts);

    ffHistEtaPtMinusAfterCuts = new TH2F("ffEtaPtMinusAfterCuts", "- Particles Vs Eta and Pt after all cuts", 200, -1, 1, 600, 0, 6);
    ffHistEtaPtMinusAfterCuts->GetXaxis()->SetTitle("#eta");
    ffHistEtaPtMinusAfterCuts->GetYaxis()->SetTitle("GeV");
    fOutput->Add(ffHistEtaPtMinusAfterCuts);

    ffHistEtaPtAfterCuts = new TH2F("ffEtaPtAfterCuts", "ch Particles Vs Eta and Pt after all cuts", 200, -1, 1, 600, 0, 6);
    ffHistEtaPtAfterCuts->GetXaxis()->SetTitle("#eta");
    ffHistEtaPtAfterCuts->GetYaxis()->SetTitle("GeV");
    fOutput->Add(ffHistEtaPtAfterCuts);

    fHistCharge = new TH1I("fCharge", "N particles vs charge", 5, -2.5, 2.5 );
    fHistCharge->GetXaxis()->SetTitle("charge");
    fHistCharge->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistCharge);

    fHistVtxX = new TH1F("fVtxX", "Vertex x coordinate distribution", 100, -0.5, 0.5 );
    fHistVtxX->GetXaxis()->SetTitle("Vx (cm)");
    fHistVtxX->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistVtxX);

    fHistVtxY = new TH1F("fVtxY", "Vertex y coordinate distribution", 100, -0.5, 0.5 );
    fHistVtxY->GetXaxis()->SetTitle("Vy (cm)");
    fHistVtxY->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistVtxY);

    fHistVtxZ = new TH1F("fVtxZ", "Vertex z coordinate distribution", 400, -20, 20 );
    fHistVtxZ->GetXaxis()->SetTitle("Vz (cm)");
    fHistVtxZ->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistVtxZ);

    fHistTPCNcls = new TH1I("fTPCNcls", "N of TPC clusters distribution", 200, 0, 200 );
    fHistTPCNcls->GetXaxis()->SetTitle("N of TPC clusters after cuts");
    fHistTPCNcls->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistTPCNcls);

    fHistITSNcls = new TH1I("fITSNcls", "N of ITS clusters distribution", 20, -10, 10 );
    fHistITSNcls->GetXaxis()->SetTitle("N of ITS clusters after cuts");
    fHistITSNcls->GetYaxis()->SetTitle("N entries");
    fOutput->Add(fHistITSNcls);

    ffHistDCAxyz = new TH2F("ffDCAxyz", "DCAxy vs DCAz", 28, -0.7, 0.7, 28, -0.7, 0.7);
    ffHistDCAxyz->GetXaxis()->SetTitle("DCA to vertex XY");
    ffHistDCAxyz->GetYaxis()->SetTitle("DCA to Vertex Z");
    fOutput->Add(ffHistDCAxyz);

    fHistChi2PerClusterTPC = new TH1F("fChi2PerClusterTPC", "Chi^2 per cluster TPC", 60, 0, 6);
    fHistChi2PerClusterTPC->GetXaxis()->SetTitle("Chi^2/Nclusters TPC");
    fOutput->Add(fHistChi2PerClusterTPC);

    fHistNch = new TH1F("fNch","N of Nch particles in event", 400, 0, 400);
    fHistNch->GetXaxis()->SetTitle("N of charged tracks");
    fHistNch->GetYaxis()->SetTitle("N of events");
    fOutput->Add(fHistNch);

    fHistNchAccepted = new TH1F("fNchAccepted","N of accepted Nch particles in event", 100, 0, 100);
    fHistNchAccepted->GetXaxis()->SetTitle("N of charged tracks accepted");
    fHistNchAccepted->GetYaxis()->SetTitle("N of events");
    fOutput->Add(fHistNchAccepted);

    ffHistEtaTPCNcls = new TH2F("ffEtaTPCNcls","N of charged particles Vs Eta and TPCNcls", 80, -2, 2, 140, 30, 170);
    ffHistEtaTPCNcls->GetXaxis()->SetTitle("#eta");
    ffHistEtaTPCNcls->GetYaxis()->SetTitle("N of TPC clusters");
    fOutput->Add(ffHistEtaTPCNcls);

    ffHistEtaPhiPlus = new TH2F("ffEtaPhiPlus", "+ particles vs EtaPhi after all cuts", 20, -1, 1, 32, 0, 2*TMath::Pi() );
    ffHistEtaPhiPlus->GetXaxis()->SetTitle("#eta");
    ffHistEtaPhiPlus->GetYaxis()->SetTitle("#phi");
    fOutput->Add(ffHistEtaPhiPlus);

    ffHistEtaPhiMinus = new TH2F("ffEtaPhiMinus", "- particles vs EtaPhi after all cuts", 20, -1, 1, 32, 0, 2*TMath::Pi() );
    ffHistEtaPhiMinus->GetXaxis()->SetTitle("#eta");
    ffHistEtaPhiMinus->GetYaxis()->SetTitle("#phi");
    fOutput->Add(ffHistEtaPhiMinus);

    ffHistEtaPhi = new TH2F("ffEtaPhi", "N particles vs EtaPhi after all cuts", 20, -1, 1, 32, 0, 2*TMath::Pi() );
    ffHistEtaPhi->GetXaxis()->SetTitle("#eta");
    ffHistEtaPhi->GetYaxis()->SetTitle("#phi");
    fOutput->Add(ffHistEtaPhi);

    ffHistEtaPhiPlus64 = new TH2F("ffEtaPhiPlus64", "+ particles vs EtaPhi after all cuts", 10, -1, 1, 8, 0, 2*TMath::Pi() );
    ffHistEtaPhiPlus64->GetXaxis()->SetTitle("#eta");
    ffHistEtaPhiPlus64->GetYaxis()->SetTitle("#phi");
    fOutput->Add(ffHistEtaPhiPlus64);

    ffHistEtaPhiMinus64 = new TH2F("ffEtaPhiMinus64", "- particles vs EtaPhi after all cuts", 10, -1, 1, 8, 0, 2*TMath::Pi() );
    ffHistEtaPhiMinus64->GetXaxis()->SetTitle("#eta");
    ffHistEtaPhiMinus64->GetYaxis()->SetTitle("#phi");
    fOutput->Add(ffHistEtaPhiMinus64);

    ffHistEtaPhi64 = new TH2F("ffEtaPhi64", "N particles vs EtaPhi after all cuts", 10, -1, 1, 8, 0, 2*TMath::Pi() );
    ffHistEtaPhi64->GetXaxis()->SetTitle("#eta");
    ffHistEtaPhi64->GetYaxis()->SetTitle("#phi");
    fOutput->Add(ffHistEtaPhi64);


    /*fHistY = new TH1F("fY", "N particles vs Y", 400, -2, 2 );
  fHistY->GetXaxis()->SetTitle("Y");
  fHistY->GetYaxis()->SetTitle("N entries");
  fOutput->Add(fHistY);*/

    /*fHistCent = new TH1F("fCent", "N events vs centrality class", 110, -10, 100);
  fHistCent->GetXaxis()->SetTitle("centrality class");
  fHistCent->GetYaxis()->SetTitle("N events");
  fOutput->Add(fHistCent);*/

    //fOutput->Add(fTrackCuts);

    //for PbPb:
    /*Int_t bins[11]={100,800,800,800,800,1600,1600,400,400,140,140};//number of bins
  Double_t min[11]={0,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-200.5,-200.5,0.195,0.195};//min values
  Double_t max[11]={100,799.5,799.5,799.5,799.5,1599.5,1599.5,199.5,199.5,1.595,1.595};//max values*/

    //for pp:       {Np_B,Np_F,Nm_B,Nm_F,Nch_B,Nch_F,Q_B, Q_F,PtMean_B,PtMean_F,Q_F*Q_B,Q_F*Q_F }

    //    TString AxisNames[12] = {"Np_B","Np_F","Nm_B","Nm_F","Nch_B","Nch_F","Q_B", "Q_F","PtMean_B","PtMean_F","Q_F*Q_B","Q_F*Q_F"};
    //    Int_t bins[12] = {  50,  50,  50,  50, 100, 100,  50,   50,  140,  140,   500,   400};//number of bins
    //    Double_t min[12]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-25.5,-25.5,0.195,0.195, -250.5, -0.5};//min values
    //    Double_t max[12]={49.5,49.5,49.5,49.5,99.5,99.5, 24.5, 24.5,1.595,1.595,  249.5, 399.5};//max values


    TString AxisNames[16] = {"Np_B","Np_F","Nm_B","Nm_F","Nch_B","Nch_F","Q_B", "Q_F","Q_F*Q_B","Q_F*Q_F", "Np_F*Np_F","Nm_F*Nm_F","Np_B*Np_F","Np_B*Nm_F","Nm_B*Np_F","Nm_B*Nm_F"};
    Int_t bins[16] = {  50,  50,  50,  50, 100, 100,  50,   50,  500,   400,   500,  500,  500,  500,  500,  500 };//number of bins
    Double_t min[16]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-25.5,-25.5,-250.5, -0.5,  -0.5, -0.5, -0.5, -0.5, -0.5, -0.5};//min values
    Double_t max[16]={49.5,49.5,49.5,49.5,99.5,99.5, 24.5, 24.5, 249.5, 399.5, 499.5,499.5,499.5,499.5,499.5,499.5};//max values


    fHistsSparse = new THnSparseF*[500];

    //here - check EtaConfiguration

        // 16, 2, 2, - 8 or 7
        // 32, 4, 2 - 15 or [12]
        // 16, 4, 2 - [7] or 5
        // 16, 6, 2 - [6] or 3
        // 16, 8, 0 - [5] or 1

    if (fNSectors == 1)
    {
        //create sparse hists only for Eta separated windows

        for (Int_t i=0;i<=fNEtaConfigs-1;i++)
            {
                TString strSparse = Form( "fEta%d", i+1);
                fHistsSparse[i] = new THnSparseF(strSparse,"Sparse hist",16,bins,min,max);

                for (Int_t d=0; d<=15;d++)
                {
                    fHistsSparse[i]->GetAxis(d)->SetNameTitle(AxisNames[d],AxisNames[d]);

                }
                fOutput->Add(fHistsSparse[i]);

            }
    }

    else
    {

        //create sparse hists for Eta-Phi separated windows

        for (Int_t j=0;j<=fNPhiSteps-1;j++)
        {
            for (Int_t i=0;i<=fNEtaConfigs-1;i++)
            {

                TString strSparse = Form( "fPhi1%dEta%d", j+1 , i);
                fHistsSparse[i+j*fNEtaConfigs] = new THnSparseF(strSparse,"Sparse hist",16,bins,min,max);

                for (Int_t d=0; d<=15;d++)
                {
                    fHistsSparse[i+j*fNEtaConfigs]->GetAxis(d)->SetNameTitle(AxisNames[d],AxisNames[d]);

                }
                fOutput->Add(fHistsSparse[i+j*fNEtaConfigs]);
            }
        }

    }

    //todo here also: if EtaConfiguration is 1,1,0 create hists for Phi separated windows only


    ffHistNpEvent = new TH2F("fNpEvent", "Nplus in one event vs EtaPhi", fNEtaBins + 2, ETA_MIN - (ETA_MAX - ETA_MIN)/fNEtaBins, ETA_MAX + (ETA_MAX - ETA_MIN)/fNEtaBins, fNPhiSteps, 0, 2*TMath::Pi() );

    ffHistNmEvent = new TH2F("fNmEvent","Nminus in one event vs EtaPhi", fNEtaBins + 2, ETA_MIN - (ETA_MAX - ETA_MIN)/fNEtaBins, ETA_MAX + (ETA_MAX - ETA_MIN)/fNEtaBins, fNPhiSteps, 0, 2*TMath::Pi());

    //fHistPtEtaEvent = new TH2F("f2PtEtaEvent", "Pt vs Eta in one event", 140, 0.195, 1.595, 16, -0.8, 0.8);

    fHistNchEvent = new TH1I("fNchEvent", "Counter of Ncharged in one event ", 1, 0.5, 1.5 );

    fHistNchAcceptedEvent = new TH1I("fNchAcceptedEvent", "Counter of Ncharged accepted in one event ", 1, 0.5, 1.5 );

    //Event level QA
    Int_t nStatBins = 8;
    fHistEventCutStats = new TH1D("fHistEventCutStats","Event statistics;;N_{events}", nStatBins,-0.5,nStatBins-0.5);
    TString gEventCutBinNames[8] = {"Total","No trigger","Centrality","No vertex","Bad vertex position","fESD not available","Mult cut","Analyzed"};

    for(Int_t i = 1; i <= nStatBins; i++)
    {
        fHistEventCutStats->GetXaxis()->SetBinLabel(i,gEventCutBinNames[i-1].Data());
    }
    fOutput->Add(fHistEventCutStats);

    //Track level QA
    fHistTrackCutStats = new TH1D("fHistTrackCutStats","Track statistics;;N_{tracks}", 7,-0.5,6.5);
    TString gTrackCutBinNames[7] = {"Total","AliESDtrackCuts","HighPt cut","LowPt cut","High Eta cut","Low Eta Cut","Analyzed"};

    for(Int_t i = 1; i <= 7; i++)
    {
        fHistTrackCutStats->GetXaxis()->SetBinLabel(i,gTrackCutBinNames[i-1].Data());
    }
    fOutput->Add(fHistTrackCutStats);

    //fOutput->Add(fTrackCuts);
}


//_______________________________________________________________________
void AliAnalysisTaskPhi::UserExec(Option_t *)

{
    // Main loop
    // Called for each event

    if (fNSectors > 0 && fNPhiSteps >0 && (fNPhiSteps > fNSectors || fNPhiSteps == fNSectors))
    {

        Float_t PhiStep = 2*TMath::Pi()/fNPhiSteps;

        Float_t PhiShift = 0.0;
        for (Int_t i=1;i<=fNPhiSteps;i++)
        {
            UserExecPhi(PhiShift);
            PhiShift += PhiStep;
        }
    }

    else
    {
        //cout<<"smth wrong with NSectors and NPhiSteps!!";
        return;
    }

}

//________________________________________________________________________
void AliAnalysisTaskPhi::UserExecPhi(Float_t PhiShift)
{


    //if ( fHaveMC )
    //	UserExecMC();
    AliESDEvent *fESD = dynamic_cast<AliESDEvent*>(InputEvent());
    if (PhiShift==0) {fHistEventCutStats->Fill(0);}     //total number of events

    if ( !fESD )
    {
        //Printf("ERROR: fESD not available\n");
        if (PhiShift==0) {fHistEventCutStats->Fill(5);}
        return;
    }

    //Trigger
    Bool_t lTrigger = kTRUE;
    if( fESD )
    {
        lTrigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
        if( !lTrigger )
        {
            if (PhiShift==0) {fHistEventCutStats->Fill(1);}
            PostData(1, fOutput);
            return;
        }
    }

    // Vertex present

    const AliESDVertex *vertex = fESD->GetPrimaryVertex();

    Bool_t lVertexPresent = ( (vertex) && (vertex->GetNContributors() > 0 ) && ( vertex->GetZRes() != 0 ) );
    if( !lVertexPresent )
    {
        if (PhiShift==0) {fHistEventCutStats->Fill(3);} //no vertex
        PostData(1, fOutput);
        return;
    }


    // Vertex in range
    Bool_t lVertexAcceptable = (TMath::Abs(vertex->GetXv()) < fVxMax) && (TMath::Abs(vertex->GetYv()) < fVyMax);
    if( lVertexAcceptable )
    {
        if( fVzMax>0 )
            lVertexAcceptable = (TMath::Abs(vertex->GetZv()) < fVzMax);
        else
            lVertexAcceptable = (TMath::Abs(vertex->GetZv()) > -fVzMax);
    }

    if(!lVertexAcceptable)
    {
        if (PhiShift==0) {fHistEventCutStats->Fill(4);} //bad vertex
        PostData(1, fOutput);
        return;
    }

    Double_t lVtxX = vertex->GetXv();
    Double_t lVtxY = vertex->GetYv();
    Double_t lVtxZ = vertex->GetZv();
    if (PhiShift==0) {
        fHistVtxX->Fill(lVtxX);//x coordinate distribution of vertex
        fHistVtxY->Fill(lVtxY);//y coordinate distribution of vertex
        fHistVtxZ->Fill(lVtxZ);//z coordinate distribution of vertex
    }

    //Printf("There are %d tracks in this event\n", fESD->GetNumberOfTracks());

    //centrality selection (only for AA)
    /*AliCentrality *centrality = fESD->GetCentrality();
   Double_t lCentrality = centrality->GetCentralityPercentile("V0M");

   // GetCentralityClass10 returns the number of 10% centrality class, an integer from 0 to 9
   // "V0M" means that we using V0 detector as the default estimator for centrality
   // other possible methods:
   // GetCentralityClass5 - number of 5% centrality class
   // GetCentralityPercentile - percentile, a float from 0 to 100 (or to the trigger efficiency)

  if (lCentrality < 0.25 || lCentrality > 80) //select centrality if needed
  {
  fHistEventCutStats->Fill(2);
  return;
  }

  //centrality classes if needs:
  Int_t lCentralityClass; //centrality class from 0-5,5-10,10-20 to 70-80

  if (lCentralityClass10 > 0 && lCentralityClass10 < 8)
  {lCentralityClass = lCentralityClass10 + 1;}
  else if (lCentralityClass10 == 0)
  {lCentralityClass = centrality->GetCentralityClass5("V0M");}*/


    // user trigger to be added (if needed)
    // user vertex cut (if needed)

    ffHistNpEvent->Reset();
    ffHistNmEvent->Reset();
    //fHistPtEtaEvent->Reset();
    fHistNchEvent->Reset();
    fHistNchAcceptedEvent->Reset();


    // ESD Track loop
    for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++)
    {
        AliESDtrack* track = fESD->GetTrack(iTracks);

        if (PhiShift==0) {fHistTrackCutStats->Fill(0);} //total number of raw tracks

        if( !fTrackCuts->AcceptTrack(track) )
        {
            if (PhiShift==0) { fHistTrackCutStats->Fill(1);} //rejectred by cuts
            continue;
        }

        if ( !track )
        {
            //printf("ERROR: Could not receive track %d\n", iTracks);
            continue;
        }


        Int_t NclustersTPC = track->GetTPCNcls();//number of TPC clusters distribution AFTER cuts
        Int_t NclustersITS = track->GetITSclusters(0);//number of ITS clusters distribution AFTER cuts
        if (PhiShift==0) {
            fHistTPCNcls->Fill(NclustersTPC);
            fHistITSNcls->Fill(NclustersITS);
        }

        // getting the track to vertex parameters
        Float_t b[2];
        Float_t bCov[3];
        track->GetImpactParameters(b,bCov);
        Float_t dcaToVertexXY = b[0];
        Float_t dcaToVertexZ = b[1];

        if (PhiShift==0) {

            ffHistDCAxyz->Fill(dcaToVertexXY,dcaToVertexZ);

            //chi2 per cluster TPC
            if (NclustersTPC!=0)
            {
                Float_t  chi2PerClusterTPC = track->GetTPCchi2()/(NclustersTPC);
                fHistChi2PerClusterTPC->Fill(chi2PerClusterTPC);
            }
        }

        Double_t lTrackPt = track->Pt();
        Double_t lTrackEta = track->Eta();
        Int_t lTrackCharge = track->Charge();

        fHistNchEvent->Fill(1); //count total number of tracks

        if (PhiShift==0) {
            if (lTrackCharge>0)
            {
                fHistPtPlus->Fill(lTrackPt);
                fHistEtaPlus->Fill(lTrackEta);
            }

            if (lTrackCharge<0)
            {
                fHistPtMinus->Fill(lTrackPt);
                fHistEtaMinus->Fill(lTrackEta);
            }

            fHistPt->Fill(lTrackPt);
            fHistEta->Fill(lTrackEta);
            ffHistEtaTPCNcls->Fill(lTrackEta, NclustersTPC);

            if ( lTrackPt > PT_MAX)
            {fHistTrackCutStats->Fill(2);}
            if ( lTrackPt < PT_MIN)
            {fHistTrackCutStats->Fill(3);}
        }

        //cut on Pt
        if( lTrackPt < PT_MIN || lTrackPt > PT_MAX )
            continue;

        if (PhiShift==0) {
            fHistPtAfterPtCut->Fill(lTrackPt);
            fHistEtaAfterPtCut->Fill(lTrackEta);

            if (lTrackCharge > 0)
            {
                fHistPtPlusAfterPtCut->Fill(lTrackPt);
                fHistEtaPlusAfterPtCut->Fill(lTrackEta);
            }

            if (lTrackCharge < 0)
            {
                fHistPtMinusAfterPtCut->Fill(lTrackPt);
                fHistEtaMinusAfterPtCut->Fill(lTrackEta);
            }


            if ( lTrackEta > ETA_MAX)
            {fHistTrackCutStats->Fill(4);}
            if ( lTrackEta < ETA_MIN)
            {fHistTrackCutStats->Fill(5);}
        }

        //cut on Eta
        if( lTrackEta < ETA_MIN || lTrackEta > ETA_MAX )
            continue;

        if (PhiShift==0) {

            if (lTrackCharge > 0)
            {
                fHistEtaPlusAfterCuts->Fill(lTrackEta);
                fHistPtPlusAfterCuts->Fill(lTrackPt);
                ffHistEtaPtPlusAfterCuts->Fill(lTrackEta, lTrackPt);
            }

            if (lTrackCharge < 0)
            {
                fHistEtaMinusAfterCuts->Fill(lTrackEta);
                fHistPtMinusAfterCuts->Fill(lTrackPt);
                ffHistEtaPtMinusAfterCuts->Fill(lTrackEta, lTrackPt);
            }

            fHistEtaAfterCuts->Fill(lTrackEta);
            fHistPtAfterCuts->Fill(lTrackPt);
            ffHistEtaPtAfterCuts->Fill(lTrackEta, lTrackPt);
        }

        //Double_t lTrackY = track->Y();
        Double_t lTrackPhi = track->Phi() + PhiShift;

        if (lTrackPhi > 2*TMath::Pi()) {lTrackPhi -= 2*TMath::Pi();}

        if (lTrackCharge > 0) {ffHistNpEvent->Fill(lTrackEta, lTrackPhi);}
        if (lTrackCharge < 0) {ffHistNmEvent->Fill(lTrackEta, lTrackPhi);}

        //fHistPtEtaEvent->Fill(lTrackPt,lTrackEta);
        fHistNchAcceptedEvent->Fill(1); //count total number of accepted tracks

        if (PhiShift==0) {
            ffHistEtaPhi->Fill( lTrackEta, lTrackPhi);
            if (lTrackCharge > 0) { ffHistEtaPhiPlus->Fill( lTrackEta, lTrackPhi);}
            if (lTrackCharge < 0) { ffHistEtaPhiMinus->Fill( lTrackEta, lTrackPhi);}

            ffHistEtaPhi64->Fill( lTrackEta, lTrackPhi);
            if (lTrackCharge > 0) { ffHistEtaPhiPlus64->Fill( lTrackEta, lTrackPhi);}
            if (lTrackCharge < 0) { ffHistEtaPhiMinus64->Fill( lTrackEta, lTrackPhi);}

            fHistCharge->Fill(lTrackCharge);
            //fHistY->Fill( lTrackY );
            fHistTrackCutStats->Fill(6);
        }
    }


    Double_t NchEvent = fHistNchEvent->Integral(1,1);
    if (NchEvent < MinMult) {
        if (PhiShift==0) {fHistEventCutStats->Fill(6);} //rejected by multiplicity
        return;
    }

    if (PhiShift==0) {fHistNch->Fill(NchEvent);}

    Double_t NchAcceptedEvent = fHistNchAcceptedEvent->Integral(1,1);
    if (PhiShift==0) {fHistNchAccepted->Fill(NchAcceptedEvent);}

    Int_t SectorSizeInBins = fNPhiSteps/fNSectors; //size of the sector in bins of Event histogram


//    //arrays of bin numbers for eta windows
//    Int_t downB[8]={9,  8, 7, 6, 5, 4, 3, 2 };
//    Int_t upB[8] = {10, 9, 8, 7, 6, 5, 4, 3 };
//    Int_t downF[8]={9, 10,11,12,13,14,15,16 };
//    Int_t upF[8] = {10,11,12,13,14,15,16,17 };


    Int_t downB[fNEtaConfigs];
    memset(downB, 0, sizeof(downB));
    Int_t upB[fNEtaConfigs];
    memset(upB, 0, sizeof(upB));
    Int_t downF[fNEtaConfigs];
    memset(downF, 0, sizeof(downF));
    Int_t upF[fNEtaConfigs];
    memset(upF, 0, sizeof(upF));


    for (Int_t n = 0; n <= fNEtaConfigs-1;n++)
    {
        if (fNSectors == 1)
        {
            // without overlapping Eta windows
            downB[n] = (fNEtaBins+2)/2 - fEtaSizeInBins + 1 - n;
            upB[n] = (fNEtaBins+2)/2 - n;
            downF[n] = (fNEtaBins+2)/2 + 1 + n;
            upF[n] = (fNEtaBins+2)/2 + n + fEtaSizeInBins;
            //cout<<downB[n]<<","<<upB[n]<<","<<downF[n]<<","<<upF[n]<<",";
        }

        else
        {
            //with overlapping Eta windows
            downB[n] = (fNEtaBins+2)/2 - (fEtaSizeInBins/2-1) - n;
            upB[n] = (fNEtaBins+2)/2 - (fEtaSizeInBins/2-1) - n + fEtaSizeInBins - 1;
            downF[n] = (fNEtaBins+2)/2 - (fEtaSizeInBins/2-1) + n;
            upF[n] = (fNEtaBins+2)/2 - (fEtaSizeInBins/2-1) + n + fEtaSizeInBins - 1;
            //cout<<downB[n]<<","<<upB[n]<<","<<downF[n]<<","<<upF[n]<<",";
        }
    }



    for (Int_t j=0;j<=fNPhiSteps-1;j++)
    {

        for (Int_t i=0;i<=fNEtaConfigs-1;i++)
        {

            if (j+SectorSizeInBins < fNPhiSteps || j+SectorSizeInBins == fNPhiSteps) {

                Double_t Np_B = ffHistNpEvent->Integral(downB[i],upB[i],1,SectorSizeInBins);//Np backward
                Double_t Np_F = ffHistNpEvent->Integral(downF[i],upF[i],j+1,j+SectorSizeInBins);//Np forward
                Double_t Nm_B = ffHistNmEvent->Integral(downB[i],upB[i],1,SectorSizeInBins);//Nm backward
                Double_t Nm_F = ffHistNmEvent->Integral(downF[i],upF[i],j+1,j+SectorSizeInBins);//Nm forward

                Double_t Nch_B = Np_B + Nm_B;
                Double_t Nch_F = Np_F + Nm_F;
                Double_t Q_B = Np_B - Nm_B;
                Double_t Q_F = Np_F - Nm_F;

                //Double_t PtMean_B = fHistPtEtaEvent->ProjectionX("name2",downB[i],upB[i])->GetMean();//Eta Window only!
                //Double_t PtMean_F = fHistPtEtaEvent->ProjectionX("name",downF[i],upF[i])->GetMean();//Eta Window only!
                //Double_t x[12]={Np_B, Np_F, Nm_B, Nm_F, Nch_B, Nch_F, Q_B, Q_F, PtMean_B, PtMean_F, Q_B*Q_F, Q_F*Q_F};
                Double_t x[16]={Np_B, Np_F, Nm_B, Nm_F, Nch_B, Nch_F, Q_B, Q_F, Q_B*Q_F, Q_F*Q_F, Np_F*Np_F, Nm_F*Nm_F, Np_B*Np_F, Np_B*Nm_F, Nm_B*Np_F, Nm_B*Nm_F};
                fHistsSparse[i+j*fNEtaConfigs]->Fill(x);

            }

            else if (fNSectors==8 && fNPhiSteps==16)  //TODO - fix this!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            {

                Double_t Np_B = ffHistNpEvent->Integral(downB[i],upB[i],1,SectorSizeInBins);//Np backward
                Double_t Np_F = ( ffHistNpEvent->Integral(downF[i],upF[i],16,16) + ffHistNpEvent->Integral(downF[i],upF[i],1,1) );//Np forward
                Double_t Nm_B = ffHistNmEvent->Integral(downB[i],upB[i],1,SectorSizeInBins);//Nm backward
                Double_t Nm_F = ( ffHistNmEvent->Integral(downF[i],upF[i],16,16) + ffHistNmEvent->Integral(downF[i],upF[i],1,1) );//Nm forward

                Double_t Nch_B = Np_B + Nm_B;
                Double_t Nch_F = Np_F + Nm_F;
                Double_t Q_B = Np_B - Nm_B;
                Double_t Q_F = Np_F - Nm_F;

                //Double_t PtMean_B = fHistPtEtaEvent->ProjectionX("name2",downB[i],upB[i])->GetMean();//Eta Window only!
                //Double_t PtMean_F = fHistPtEtaEvent->ProjectionX("name",downF[i],upF[i])->GetMean();//Eta Window only!
                //Double_t x[12]={Np_B, Np_F, Nm_B, Nm_F, Nch_B, Nch_F, Q_B, Q_F, PtMean_B, PtMean_F, Q_B*Q_F, Q_F*Q_F};
                Double_t x[16]={Np_B, Np_F, Nm_B, Nm_F, Nch_B, Nch_F, Q_B, Q_F, Q_B*Q_F, Q_F*Q_F, Np_F*Np_F, Nm_F*Nm_F, Np_B*Np_F, Np_B*Nm_F, Nm_B*Np_F, Nm_B*Nm_F};
                fHistsSparse[i+j*fNEtaConfigs]->Fill(x);

            }

            else {
                //cout<<"smth wrong with eta-phi configurations!";
                return;
            }


        }
    }



    //fHistCent->Fill(lCentrality);
    if (PhiShift==0) {
        fHistEventCutStats->Fill(7);
    } //analyzed

    PostData(1, fOutput);

}



//________________________________________________________________________
void AliAnalysisTaskPhi::Terminate(Option_t *)
{
    // Draw result to the screen
    // Called once at the end of the query

    fOutput = dynamic_cast<TList*> (GetOutputData(1));
    if (!fOutput) {
        //printf("ERROR: Output list not available\n");
        return;
    }

}
